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ABSTRACT 


The object of this thesis 1s to examine, vila Simulation, 
the properties of a class of single-server, first-come 
first-served queues in which the service times and inter- 
arrival times, while still marginally exponential, are 
auto-~correlated and cross-correlated. The correlation is 
introduced by letting the service and interarrival sequences 
be EARMA-type processes, where EARMA stands for exponential 
autoregressive, moving average sequence. An extension of 
these ideas brings in the cross-correlation. The waiting 
times in the dependent queue are compared to the waiting 
times (with known distribution) in the independent M/M/1 
queue uSing data analytic and formal statistical methods. 
Varlance reduction techniques are also studied; these use 
distributionally known aspects of the M/M/1 queue (Simulated 
with the same exponential sequences as the dependent queue) 


to control the unknown aspects of the dependent queue. 
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ioe sNTRODUCTION 


The object of this thesis is to examine, via simulation, 
the properties of a class of single-server first-come first- 
served queues in which the service times and interarrival 
times, while still marginally exponential, are auto-correlated 
and cross-correlated. The correlation is obtained by making 
the service and interarrival time sequences multivariate EARMA 
processes. An ancillary aspect of the study is the investi- 
gation of the efficiency of a number of control variable 
variance reduction schemes. These schemes are made possible 
by the simplicity of the structure of the multivariate EARMA 
processes, and the fact that the null case (no correlation) 
gives the M/M/1 queue, whose properties are known analytically. 

An immediate problem which arises in the simulation is 
that it 1s not known how far out along the sample path the 
Simulation must be carried in order for the steady state 
properties of the queue to hold. This 1s a general problem 
which arises in queuing Simulations and, in order to cope with 
it, the simulations are run in 500 parallel replications. 

This allows one to use many methods of time series analysis 
and data analysis to obtain graphical verification in a se- 
quential manner of the approach to a steady state. The repli- 
cations also allow for examination of the complete steady 
State waiting time distribution in the queue, and not merely 


the mean waiting time. 





ale 


While all of these aspects of the queuing Simulation 


interwoven, we can separate them out for purposes of 


discussion as follows. 


(@) 


ee) 


An elementary aspect of the queue which can be examined 
is the mean waiting time. The most efficient (smallest 
variance) estimate of this is the sample path average 
which is again averaged over replications. This repli- 
cation of course reduces the variance by a factor of 
m, where m is the number of replications. Furthermore 
the existence of replications allows one to look at the 
complete distribution of the sample path average esti- 
mate. This distribution should be normally distributed, 
possibly even before a steady state has been reached. 

The object of this part of the study is to see how 
the mean steady state waiting time in the queue varies 
with the traffic intensity, t, and the correlation 
parameters, and in particular, how this mean waiting 
time differs from the waiting time in the M/M/1l queue 
(null case). 
A more detailed aspect of the queue which has to be 
examined is the complete distribution of the steady 
state waiting time. In the M/M/1l queue this is 
exponential, given that the waiting time is greater 
than zero. 

The first problem here is to determine that one 
1s actually looking at the steady state waiting time. 


This is done by examining the empirical distribution 





fr. ) 


function of the waiting times across replications 


at times ny and na (n. > ny 


refers to the customer number. Of course the two 


), where the index n 


samples from "ny and fe each of size m, will be 
correlated and this makes standard statistical methods 
for comparing samples invalid. However the correla- 
Eien between these samples can be calculated in order 


to determine that n. is sufficiently large relative to 


2 


n, for the correlation to be negligible. In actual 


1 
fact the comparison is done by looking at box plots 
of the ws at successive values n, <n, < 7, Se epee 
This method iS graphic and allows the simulator to 
see how the steady state 1s being approached. 

A third problem which is involved in looking at the 
complete distribution of the waiting time is that a 
sample of size m = 500 is not sufficient to really 
get a good estimate of the steady state distribution 
of the waiting time. Thus successive samples along 
the parallel sample paths are combined, where it is 
determined that the samples are far enough apart to 
be approximately independent. 

Again the object here is to see how the complete 
distribution of the steady state waiting time differs 
from the exponential distribution (null M/M/1 case) 
as traffic intensity, t, and correlation parameters 


vary. In particular heavy traffic theory says that 


this distribution should be approximately exponential 


10 





as the traffic intensity parameter approaches l. 

Tae aifrrcul ty nin weaedzptine this, and the heavy 
traffic approximation to the mean waiting time, is 
that the time to reach the steady state becomes very 
long aS t > l. An additional question that arises 
then is whether it 1s possible to start the simulation 
in such a way (approximate steady state conditions) 

SO as to Speed up this convergence. It should be 
remembered, that the queue with correlation is not 
regenerative and no initial ponacen tone for stationarity 
are known or are likely to be found. 

Variance reduction is always helpful in a simulation 
of this kind since real detail is always masked by 
Sampling fluctuations l1.e., the variability in the 
estimate of the unknown mean waiting time. Because 

of the simple probabilistic-linear structure of the 
EARMA processes, running on M/M/1 queue in parallel 
with the correlated queue uSing common exponential 
variates provides a very good control variable for the 
Waiting time if the correlation, in a rough sense, 

in the EARMA queue, iS small. However, as this 
correlation increases, the amount of control decreases. 
Thus several other control variables are examined and 
combined into a multiple control for the correlated 
queue. In this way it is hoped to obtain variance 


reduction for any values of the correlation parameters 


ae 





in the EARMA queue. The actual degree of control 
Vitchmcctimeemctttalncdets Cxatined, using the repli- 


cated simulation, as a function of the parameters 


of the queue. 


IL 





II. THE EXPONENTIAL AUTOREGRESSIVE-MOVING 
AVERAGE PROCESS, EARMA(1,1) 

The first order moving average process and the first- 
order autoregressive processes were combined by Jacobs and 
Lewis (1977) to form the EARMA(1,1) Sequence. We first 
describe these two first order sequences, and then the method 


of combining them. 


A. THE EXPONENTIAL AUTOREGRESSIVE PROCESS, EAR(1) 
The standard linear, first-order autoregressive model for 
a stationary sequence of random variables {X,} is defined 


by the equation 


cra .1) Re = lO Kis cme Oe tf) Ot 2 key 


where 9 iS a constant which is less than 1 in absolute value 
and meat 1s a sequence of independent and identically dis- 
tributed random variables. Gaver and Lewis (1978) showed 

that if the {X,} sequence were to have an exponential marginal 
distribution with parameter >’, then the parameter o0 should be 
greater than or equal to zero and less than one, and E. 

should be zero with probability o0 or an exponential (A) 


random variables, Ea, Vite cOndobeiiy lo. Thus 


becomes 


3 





wlth probability 0 


(oer. 2) X 


OX. + E with probability l-o, 


me O, +1, +2, ... where {E. is an identical independent 
distribution (1.1i.d.) Sequence of exponential (A) random 
variables. Note that for this EAR(1) model the distribution 
of the 5 depend on o, the multiplicative weight of Xela: 
Also or is not an absolutely continuous random variable. 
Thus standard results (Mallows, 1968) to the effect that the 
X; in an autoregressive process become approximately normal 


as 9 > 1 do not hold. In fact in the EAR(1) process (II.2), 


the X, eaeme geonenclabiy distributed for 0 < p < l. 


B. THE EXPONENTIAL MOVING AVERAGE PROCESS, EMA(1) 
The first order moving average exponential process EMA(1) 


was developed by Lawrance and Lewis (1977) in the form of 


oe. 3) X 


Sa GP The W.ep. (1-8), 


1 = 0, +1, +2, ... , where 8 is greater than or equal to zero 
and less than or equal to one, and iE; 4 is again a sequence 
of 1.1.d. exponential (i) random variables. (This is the 
backward case — a forward case 1s also possible.) The X's 


have an exponential marginal distribution and are only serially 
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dependent for lag one; this model is highly tractable and 
a full account of the statistically useful properties was 
obtained by Lawrance and Lewis (1977). The forward model 


combines E. with E. instead of with E. 
al it i- 


1 1° 


C. THE EXPONENTIAL AUTOREGRESSIVE-MOVING AVERAGE 
PROCESSES EARMA(1,1) 
The EAR(1) model of Section IIA and the EMA(1) model 
of Section IIB are combined in the following form to give 
an EARMA(1,1) process. This is a non-Markovian model which 
contains the EMA(1) and EAR(1) models aS a Special case. The 


defining equations for the EARMA(1,1) process are 


BE, W.P. By 
eum. 4 ) x; = 
SE. + Asya Weer 6, 
(Cee ol = 0; tl, +2, ) 
where 
pA; _» WaDie 05 
ey 
pA; _» SF ES 1 5 joa ler 


and {E, } is, aS usual, a sequence of i.i.d. exponential (A) 
random variables. Thus instead of E. being combined with the 


exponential random variable Esq it is combined with an 


Ls 





exponentially distributed random variable which is an EAR(1) 
combination of Esq! E,i91 Eyiagr cee The Sena correla- 
tions for the EARMA(1,1) process are given in Jacobs and 


Lewis (1977) as 


mmeeec(o, 8) = 8(1-8) (l-p) + (1-8). 


Note that when o = 0 the EARMA(1,1) process reduces to an 
EMA(1) process; if 8 = 0 it is the autoregressive EAR(1) 
process, and when 8 = 1 it is the usual Sequence of 
exponential (A) random variables. We will say that Xs 


is autoregressive over Ew Ely: 
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MVtton OC Max D STRUCTURES EARMA(T,T) 
IN MODELLING QUEUES 
Consider for simplicity a queue with a single input 
stream and a Single server and a first-come-first-served 
(FIFO) service discipline. Let Sir 1 Saal, opt ee Y 
emote the Service time for the ith arrival, and let Xs, 
Mol, 2, --. , denote times between arrival of the ith and 
(i-l)th customers (interarrival times). As iS uSual we aSsume 


that the first customer (with service time S.) arrives at 


0 
time zero and finds the queue empty. 

If the Sa hi and {X,} Sequences are i.i.d. exponential 
random variables with parameters \» and a respectively, we 


have M/M/1 queue. 


Now let 


E; pee inencemexponenecials (i); 1= 0, 1, 2, ««s , 


Es Demi saw nexponenuaal (a); -1-= 0, 1, 2, 


We want to model queues with correlated (autocorrelated 
and/or cross-correlated) service and inter-arrival times, 
the service and inter-arrival times both having marginally 


exponential distribution (MDxMD/1 queue). Exponential 


marginals iS a common assumption and with it the MDxMD/1 


queuing model includes the M/M/1l queue as a Special case. 


JED 





The dependence in the arrival and service processes 
is created here by defining S, and X, asa bivariate depen- 
dent sequence of random variables with exponential marginal 


distributions. This is done by letting {S,} be EARMA(1,1) 


over E,, (afA)es, (A/A) Es yy 5 Sp) these 
So = Ey 
BE, w.p. 8B; A, = (a/r) ey; 
Sy = 
BE, = Ayr Wepre (loeb). 
SE. Wis. BF pA, Wo ae 
S5 = A, = 
O 
3E. “ Aoy wisp (hos); pA, + TE. W-P- (1-0); 
BE., DEB owe 1s) : 
al oe aL Wie O:; 
S . = Ne pee 
ae aL 
BE. + A., w.p. (1-8); a rs 
a als PA.) a= Te, W-P- (1l—p) ; 
We also let 
mA = €.3 Leas 2, 3; ; 
L ab 
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although one could further make the X's correlated. With 
the present assumptions the input process is still a Poisson 
process, aS in the M/M/1 queue. The cross-coupling between 
the sequence {S.} and {X, } is apparent from the structure. 
Note that the S.'s are now correlated because of the cross- 
Seupling. The iS, 1s a sequence of dependent exponential 
random variables, but not an EARMA(1,1) sequence. It is a 
type of pseudo-EARMA(1,1) sequence. 
The interpretation of this queuing model is that the 
server tends to speed up if the queue gets long (in the 
past). Of course he also Slows down when the queue gets 
Short. In the case where 9 = 0 and g = 0 then the correlation 
between {S5} and iX,} will be equal to one, i.e., o 
S; = (a/A)X, Also when 8 = 1, and for any value of o, 
the queuing model will be identical to the M/M/1 queue. Thus 
the M/M/1 is included as a special case. Some more complicated 
schemes for coupling the service and inter-arrival processes 
are given in Lewis and Shedler (1978). Jacobs (1978a, 1978b) 
has studied heavy traffic approximations for EARMA-type queues. 
Other attempts to model M/M/1 (FIFO) type queues with 
correlated service and interarrival times have been made. 
None of them, however, have given fixed (exponential) 
marginal distributions for S; and Xsy or cross-coupling between 
these sequences. There are, of course, Markovian schemes 
based on birth-death representations which give FIFO queues 


in which the service times and interarrival times are correlated 


ike, 





(Cox and Smith, 1961), but the scheme is quite different 
from that given here. In particular the interarrival and 


service times are not exponentially distributed, and this 


may not conform to what is observed in practice. 
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IV. THE SIMULATION MODEL 


A. INTRODUCTION AND NOTATION 

Since both queues we are considering, the M/M/1l queue 
and the MDxMD/1, are single-server, first-in-first out queues, 
the waiting times Wa im SOeweare given by the recursion 


equation 


W = max(W +S —- Xx Oe es n= 0 4.2, tice 
n n 


n+l’ 
These equations are used to generate successive wis is ene 
Simulation. This is the only aspect of the queue which will 
be considered. To have looked at quantities like the number 
of customers in the queue at departure times would have com- 
plicated the programming and taken away from the primary aim 
of the thesis, which is to explore the effect of correlation 
on the queue and to use some relatively modern statistical 
methodology to do this. 

We denote the waiting time of the nth customer to 
arrive in the jth realization of the queue by W,(J)- The 
usual sample path estimate for the mean, E(w), of the 


limiting distribution of waiting times, BUX), where 


Fo(x) = lim Piw. < xi , 
W Beer n 


Za 





is the sample path average 


n 

a 

RSS Ss ) a) 
L=0 


which can be shown to converge to the mean E(W) as n > ~. 

The M independent sample path realizations We) 
emi, 2, ---, } = 1,2, .---, M can be used to obtain estimates 
of 

Meee the distribution of the estimate wi (3)- This should 

be normal for large n. This can be examined from 
the sample Wi (3), o— 1 eee. «olimeaddition the 


mean of this sample 


= M M n 
—- 4 ee Ee . 
neo) eee) |) Ww, (3) 
j=1 j=l 2=0 


is a grand mean estimating E(w) whose variance is 
computable from the sample ye i ee ssp 
Thus 


Var (W = = (sample variance of W, (3) 's) 


g? 


M 
a ree | = 
=) Cres ) (Wij) - Wi”) 
j=1 
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Cree the -rstcrabution of the (not only itS mean) from the 
sample Soe ia, ..., M. 9 )Since this 1s a random 
sample, all the classical methods are available, 
Seno eistog@als,;enormal pilots, empirical ¢c.d.f's. 

(3) The correlations between successive waiting times, 
such as a): De ees Sees 1S 


a: Z 
given as 


divided by the sample standard deviations of the two 


samples, where 


= 

| 
=| bh 
ee 
= 

‘a 


1S an average across replications. 


Pee Pe ROGRAM STRUCTURE 
l. General 
To simulate this particular EARMA(1,1) with cross- 
correlated service time and inter-arrival time, a FORTRAN 
program has been written to calculate the mean waiting time 
W, fd) - It also calculates the mean of various statistics 
for the M/M/1 queue uSing the same exponential deviates that 


are used to generate the correlated queue. These statistics 


Z3 





are used to control the estimates from the MDxMD/1 queue for 
purposes of variance reduction. The maximum number of 
customers which the program can handle for each replication 
is n = 10,000. This computation is divided into a maximum of 
10 steps and the samples of waiting times at these points are 
stored in a file. Thus for example, if we use the maximum 
size n = 10,000, we have ten data sets W 


1000 '3)* “2999 3? 
a Wi9,000'3)* ia ee.) Tous the time evolution of 
the queue can be studied. Moreover there is a facility to 
restart the simulation to get Wi1,000'4)" Wo 000 19) Bees 

This way the evolution of the queuing process can be studied 

as far out as needed, e.g., to n = 270,000, or n = 1,000,000 etc. 
Another program will take care of reading from the output file 
to do data analysis, 1.e., box ,plot of waiting times at 

various steps, average waiting times a) and their statistics. 


Zee MAINI Program 


a) The variables are to be read in as follows: 


N = Number of arrivals to be generated 
(max 1s 10,000) 


RS, RS = Arrival, Service rate. 

BETA, RHOS = 6, 0 

KN = Number of the run 

KS = 1 irmdinast wie sOumcenaton = 0,1,....,N 
0 Otherwise (1.e., restart) 

NREP = Number of replications (max 500) 

NSTEP = Number of steps in this run 


(maximum of 10) 
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NPAST 


1h IES: 


IBRTA, IRHOS 


= Number of arrivals of last run 
(-ee last Fun Stopped at n = 120,000) 


Seeds to generate exponential random 
variable of Arrival, Service times 


Seeds to generate uniform random 
variables to generate probabilities 
in the EARM(1,1) process of moving 
average and auto-regresSive. 


b) The output will be stored in two separated files. 


a) 


For continuation of the generation of the waiting 


times of each replication to the next run (1.e., the restart), 


the variables stored are as follow: 


ARIV 


SERV 


SNEW 


WAIT 


WBAR 


ALAST 


ae, 


TR 


the correlated 


Il 


Correlated queue 


cumulative inter-arrival time (these are 
the same for both queues) 


eumulative service time 


Soniiee ft imeuoretenme last customer in this 
sanbbel 


waiting time of the last customer in this 
run 


cumulative waiting time 


Auto~regressive component of the inter- 
arrival time of the last customer in this 
run 


Number of customers who have arrived since 
the last zero waiting time, 1.e., backward 
recurrence time in the renewal process of 
times of emptiness in the queue 


Number of zero waiting times up to the 
Corminaelonsoornme Of this run. 


Uncorrelated queue 


The variables are defined as the same as for 


queue; in the program the quantities are 


2S) 





suffixed by M, e.g. WAITM is the waiting time for the 
uncorrelated queue and WAIT the waiting time for the 
correlated queue. 
ii) For analyses and displays of the data in 
each run and each step the data, including the data with 
control variables added, are stored in another file in the 
sequence of waiting time in each replication of the corre- 
lated queue, average waiting time in each replication of 
the correlated queue, waiting time in each replication of the 
uncorrelated queue, average waiting time in each replication 
of M/M/1l queue, waiting time (controlled) and average waiting 
time (controlled) of the correlated queue. 
3. Subroutine CARMA 
This subroutine is used to create a string of service 
times, Sar which are cross-correlated to the inter-arrival 
time Kos and also to create the moving average structure. 
4. Subroutine CONVAR, COPYSM, COPYMV and MPROCV 
These four subroutines are used to compute the 
Statistics from the M/M/1l queue which are used to control 
Bye data from MDxMD/1 queue for purposes of variance reduction. 
5. MAIN2 Program with Subroutine COMPAR, FILL and STAT 
These programs and subroutines are used to display 
and analyze the data for each run. They produce the figures 


given in this thesis. (See Appendix II.) 
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a. Description of Box Plot 

Subroutine COMPAR will produce vertical box 
plots (McNeil 1977) parallel to each other and write out 
the minimum and the maximum value of all the data in the 
array. 

Each box plot is obtained by first calculating 
the lower and upper quartiles and the median of the batch 
of numbers. 

Where the median of a batch of numbers is the 
value for which half of the numbers in the batch are larger 
meaeait are Smaller, the lower quartile is the value that 
divides the batch into two parts, with 1/4 of the numbers 
below this value, and 3/4 above it. Similarly the upper 
quartile is the value for which 3/4 of the numbers are 
below it and 1/4 above it. 

The narrow rectangular box with ends corres- 
ponding to the lower and upper quartiles contains the median 
point (*). The height of the box, the interquartile 
distance, D , 1s then measured out on each side of the 
quartiles. Then the lowest and highest data points that 
fall between these quantities are also marked by crosses (X). 
Finally, any numbers whose positions are outside these crosses 
are marked with circles (0), those more than 1.5 inter- 
quartile distances outside each quartile getting the 
symbol (@). The digit besides the "outlier" positions is_ 
to indicate the number of data points that are too close 


to be marked separately. (See Figure l .) 
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V. ANALYSIS OF SIMULATION RESULTS 


me |6CLNTRODUCTION 

The object of the Simulation was to examine the effect 
of the cross correlation, indexed by the coefficients 0 and 
8B in the EARMA(1,1) process, on the stationary waiting time 
W for various values of the traffic intensity = Both the 
mean and distribution of W were compared to the known mean 
and distribution of the waiting time in an M/M/1l queue 
(Same t and same mean service time). 

To perform this comparison we initially simulated this 
meael for values .25, .50, .75, .90 of the taffic intensity 
memecor Values 0.0, .25, .50, .75, & .90 of the autoregressive 
parameter (0), and the single value 0.50 for the moving aver- 
age coefficient 8) for each combination of t and 0. These 
values were chosen to limit the amount of computing to a 
physically feasible range, with the idea that higher values 
of 90 and t would be examined if time permitted. The number 
of customers, n, until steady state is reached goes up drama- 
tically as 9 and t increase. The Simulations are run in 500 
parallel replications, and the 10,000 arrivals for each run 
are divided into 10 steps of 1,000 arrivals. 

Successive groups of 10,000 arrivals were run for each 
combination of 0 and t until it was determined from an analy- 
Sis of the output that the simulation had reached a steady 


State. This determination was made from box-plots of samples 


2g 





en’: Ma nes ele, My, Lor no 


a and Woe and from formal statistical comparisons of the 


> Nyy Peon DOS OF 


two samples with n., large enough compared to n, that the 


Z iE 


correlation between the two samples is negligible. 


Bee GoLlLIMATION OF TIME TO STEADY STATE 

For each run, the distribution of the sample path average 
estimate W_ (3) for each step is observed by the method of 
box plots. By comparing these plots along the n axis at 
steps of arrivals of 1,000, convergence to steady-state and 
normality cei be determined. Figures la-g show the conver- 
gence to steady state of the M/M/1 queue at t = .75. For 
the EARMA-type MDxMD/1 GUeWCEWretee = .75, p = .75, 6 = .50, 
the box plots are shown in Fig. 2a-g. The number of arrivals 
at the estimated steady Sete. rOm@muene M/M/l queue and. for 
the correlated queue for any combination of pop and t are 
shown in table l. 

In Figures 1 and 2 note that the variance of ah is 
decreasing as l/n, so that there is a continual shrinkage 
Peome.ert to right in the plots. The symmetry and lack of 
extreme values gives an indication of normality of the esti- 
mates. Note that the scales on successive figures change 
through la to lg etc. This 1s because the plots are not 
commensurate. Range of values is indicated by the values 
Geesop and bottom on the left, e.g., 2.3460522 and 0.3639930 


in Figure la. 
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Correlation, 


tethrre = tntcensity, tt 


40,000 Se O10 


50,000 50,000 


40,000 50,000 50,000 70,000 


40,000 50,000 10 8010 100,000 


60,000 70,000 100,000 OOy Coe 


Table l. 





Number of arrivals that the simulation 
run indicates are needed approximately 

to reach steady state (for §=0.5). These 
are read off of plots such as those given 
aoe —Oey 5 eo = NOe7> and 6B = 0.5 in 
Figures 2a-g. 
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Figure la. Box plots derived from 500 replications, j=1,...,500 of the 


sample path estimate W,(j) of the mean waiting time. The 
mean of these,W,, is an overall estimate of the mean waiting 
time. The medians of the Wy(j) are given by the * in the 
box plot. 


Uncorrelated queue (M/M/1) 


N 
t 


1,000 to 10,000 in steps of 1,000 
0.75, 1/E(S) = Service rate = 4.0 


2 


RIA) 
AkARY NO 


>< OIE 





1.€272485 


a. 


A) 


A) 











a a 
a a 
a) a 
@ r} 
a 
a ry 
a 
0 a3 7) rs) a a 
Q 0 a2 a @ 
C3 Oe 0 ' a 
04 O 93 Q3 Ore a2 a3 
C2 0 03 02 Or; GZ CZ C 
Oz Q3 02 9 J3 92 2) 
xX 02 O G3 C C4 Ge a4 
X x X X 05 Oz 06 C4 
xX x Oe Gs 
ff 
| | ) 
-+- -+— | | 
| 1] [| | | [| ~+— | [| 
; * } = ~ ste 4 x + 
! | 
| | 
| | al 1 | | 0 =4- 
-+- —+- -+- -+- -+- -+- -+- | 
~+- | | : 
| 
! , | 
| | x X x 
| X x Q 32 C2 
| X X X 0 0 
x 03 02 02 0 Oz 0 C 
C2 92 8) 2 Q ) Nn 
0 ) 0 C 0 C 
: 0 a a 
Oe ) aq a qa) a 
a 
a a 
O27 CC 161 


Figure 1b. 


Box plots derived from 500 replications, j=]1,...,500 of the 
sample path estimate Wy(j) of the mean waiting time. The 
mean of these, W,, is an overall estimate of the mean waiting 


time. The medians of the Wn(j) are given by the * in the 
box plot. 


Uncorrelated queue (M/M/1) 


N = 11,000 to 20,000 in steps of 1,000 
t = 0. is 1/E(S) = Service rate = 4 0 


i] 


38 


ee ee | 


@) 


he ne ee 5 OWI ODOR) 
mini hd 


+ 


ve 
pee 


HYAyOoo aOrn6§ 
Nh 





0.8985164 
a 


Bap 


xODOOOOADMO A) 
b SNMP EN) 


(| 
Sg 
j ee ee ( 





a 
0 
0 03 a q@ a 
05 02 0 Q3 a 
0 a) Q2 ¢ 8) Q 
03 04 O02 C2 C4 3 C C2 
95 96 04 02 0 O2 Q3 0 
5 04 O5 oe Os Cz 0 Q3 C 
X 04 G3 C2 os Q2 O7 
X x X X X x 0 
| | C5 
Xx 
| 
i 
| 
+ + 
+ 


+ 


| | 


at 
i———— 1 
% 
Ua | 
+ 
{—_—— «I 
+ 


j ge ge lta ae | 
+t 
| a eel ] 


i— 








+ aa + ~ | + 
: | | X 
x xX x X x C 
X x X x Ole! Q3 C2 0 Q Oz 
a5 9 33 32 92 02 4) 9 C 
| D2 6 
0 0 
C2 9 92 0 C 02 a 
02 0 9 C a 
e a) al 
a a) 
a 
a 
ad 
az 
Heong =t7S 
Figure lc. Box plots derived from 500 replications, j=l1,...,500 of the 


sample path estimate Wn(j) of the mean waiting time. The 
mean of these, W,, iS an overall estimate of the mean 
waiting time. The medians of the Wy(j) are given by the * 
in the box plot. 

Uncorrelated queue (M/M/1) 


N = 21,000 to 30,000 in steps of 1,000 


t = 0.75, 1/E(S) = Service rate = 4.0 
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Figure ld. 


Box plots derived from 500 replications, j=1,...,500 of the 
sample path estimate W,(j) of the mean waiting time. The 
mean of these, Wn, is an overall estimate of the mean 
waiting time. The medians of the W,(j) are given by the * 
in the box plot. 
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Figure le. Box plots derived from 500 replications, j=1,...,500 of the 


Sample path estimate W,(j) of the mean waiting time. The 
mean of these, Wy, is an overall estimate of the mean 
waiting time. The medians of the W,(j) are given by the * 
in the box plot. 

Uncorrelated queue (M/M/1) 

N = 41,000 to 50,000 in steps of 1,000 

0.75, 1/E(S) = Service rate = 4.0 
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Figure lf. Box plots derived from 500 replications, j=1,...,500 of the 


sample path estimate Wy(j) of the mean waiting time. The 
mean of these, W,, 1s an overall estimate of the mean 
waiting time. The medians of the W,(j) are given by the * 
in the box plot. 
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Figure 1g. Box plots derived from 500 replications, j=1,...,500 of the 
Sample path estimate W,(j) of the mean waiting time. The 
mean of these, Wy, iS an overall estimate of the mean 
waiting time. The medians of the W,(j) are given by the * 
in the box plot. 

Uncorrelated queue (M/M/1) 
N = 61,000 to 70,000 in steps of 1,000 
t = 0.75, I1/E(S) = Service rate = 4.0 
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Figure 2a. Box plots derived from 500 replications, j=1,...,500 of the 


sample path estimate W,(j) of the mean waiting time. The 
mean of these, Wy, iS an overall estimate of the mean 
waiting time. The medians of the Wy(j) are given by the * 
in the box plots. 


Correlated queue (MDXMD/1) 
N = 1,000 to 10,000 in steps of 1,000 


te-70n7on f= 0.5, 0 = 0.7/5. Service rate = 4.0 
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Figure 2b. Box plots derived from 500 replications, j=1,...,500 of the 
sample path estimate W,(j) of the mean waiting time. The 
mean of these, W,, is an overall estimate of the mean 
waiting time. The medians of the Wy(j) are given by the * 
in the box plots. 
Correlated queue (MDxMD/1) 

11,000 to 20,000 in steps of 1,000 

0.75, 8 = 0.5, 9 = 0.75. Service rate = 4.0 
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migume 2c. 


Box plots derived from 500 replications, j=1,...,500 of the 
sample path estimate Wn(j) of the mean waiting time. The 
mean of these, Wy, is an overall estimate of the mean 
waiting time. The medians of the Wy(j) are given by the * 
in the box plots. 

Correlated queue (MDxMD/1) 

N = 21,000 to 30,000 in steps of 1,000 

t = 0.75, 8 = 0.5, p = 0.75. Service rate = 4.0 
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Figure 2d. Box plots derived from 500 replications, j=1,...,500 of the 
Sample path estimate W,(j) of the mean waiting time. The 
mean of these, W,, iS an overall] estimate of the mean 
waiting time. The medians of the Wy(j) are given by the * 
in the box plots. 

Correlated gueue (MDXMD/1) 
N = 31,000 to 40,000 in steps of 1,000 
t = 0.75, 8 = 0.5, p = 0.75. Service rate = 4.0 
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Figure 2e. Box plots derived from 500 replications, j=1,...,500 of the 
Sample path estimate Wn(j) of the mean waiting time. The 
mean of these, Wy, 1S an overall estimate of the mean 
waiting time. The medians of the Wn(j) are given by the * 
in the box plots. 

Correlated queue (MDxMD/1 ) 
N = 41,000 to 50,000 in steps of 1,000 
t = 0.75, 8 = 0.5, p = 0.75. Service rate = 4.0 
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Figure 2f. Box plots derived from 500 replications, j=1,...,500 of the 


sample path estimate Wn(j) of the mean waiting time. The 
mean of these, Wy, is an overall estimate of the mean 
waiting time. The medians of the Wn(j) are given by the * 
in the box plots. 


Correlated queue (MCXMD/1) 


N 
2 


51,000 to 60,000 in steps of 1,000 
0.75, 8 = 0.5, 9 = 0.75. Service rate = 4.0 
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Figure 2g. Box plots derived from 500 replications, j=1,...,900 of the 
Sample path estimate Wy(j) of the mean waiting time. The 
mean of these, Wn, iS an overall estimate of the mean 
waiting time. The medians of the W,(j) are given by the * 
in the box plots. 

Correlated queue (MDXMD/1) 
N = 61,000 to 70,000 in steps of 1,000 
t = 0.75, 8 = 0.5, 9 = 0.75. Service rate = 4.0 
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C. ESTIMATED MEAN WAITING TIMES 

To compare the change of the estimate steady state mean 
waiting time from the M/M/1 queue as po is changed, Table 2a 
shows the value of mean waiting time, a and the estimate 
Standard deviation of those means. The mean service rate 
is fixed at 4.0. 

The mean waiting times of the correlated queue for any 
fixed value of 0 and t, are plotted in Figures 3a and 3b. 
respectively. These show that the mean waiting time increases 
dramatically as t > l, just as for the M/M/1 queue (8 = 1.0), 
for which the mean meres cine increases as 1l/(l-t). A 
Similar drastic increase in the mean waiting time because of 
the increase in cross-correlation, measured by Deis (apparvene 
in Figure 3b, when t is large. Note that when o is close 
to one there is very long term correlation in the queuing 
process; however, 0 = 1 is not allowed because the system is 
not ergodic. 

The ratio of the estimated mean waiting time, WT £Or 
the correlated queue to that for the M/M/1l queue which are 
tabulated in Table 2b, are also plotted in Figure 4a and 
Figure 4b. In particular, Figure 4a shows that the correla- 
tion causes a decrease in mean waiting time for o small, 
and an increase when 0 is large. The latter effect probably 
occurs because the service times become highly autocorrelated 
as 9 increases. 


The shape of the plots in Figure 4a suggested the possi- 


bility of fitting the mean waiting time to a function of 
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o0/(1-e) and t/(l-t) so as to be able to predict E(W) for 
large values of t. The form of the independent variables 
were. suggested by the fact that a log transform showed 
MmE(W) to go up approximately as In(l-o) for large p, and 
the fact that for small 0, E(W) is smaller than E(W) in 
the M/M/1 case. Also the form t/(l-t) was suggested by the 
M/M/1 theory. 

Only one value of 8 was used in the Simulations, 8 = 0.5, 
so that the results did not reflect this variable. The rough 
fit obtained for t close to one was, for the MDxMD/1 queue 


with parameters t,p,E(S). 





7.1) ain (en). ee Sync pe ee 


1-t Tn aay 


One use for this formula would be to predict E*(W) for 
high t so that a simulation could be started with Wo exponen- 
tial with mean E(W_) (the Duebabtlity Of a Zero Waiting time 
can be neglected for high t). This reduces the transient 
effect, as we will see later. 

meemuta (V.1) predicts that for t = 0.99, pop = 0.25 


and E(S) = 0.25, the mean waiting time would be 


Boe enae, tf = 0.99, 6 = 0.25, E(S) = 0.25 


The mean for the M/M/1 queue is given as 
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A simulation was performed for this traffic intensity 
mer— .99), starting with Wo = 0; the simulation was found to 
converge to steady-state for n = 280,000 and then (see 


Begures 12a and 12b) 


W280,000 oe Pa Od 


with 


= | 


Serle eee i O1e, 


280,000 


This is in reasonable agreement. 
The effect of using this value for E(W)) in the simulation 


is discussed in Section V.D.. 
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Ratio of the Estimated mean W, of 
the average waiting time in the 
correlated queue to the known 
average waiting time of the 
uncorrelated (M/M/1) queue (8 = 0.5) 
for various values of o0 and t. 
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Figure 4a. Smoothed, estimated mean waiting time in the correlated 
MDxXMD/1 queue divided by the known mean waiting time in 
the M/M/1 queue as a function of traffic intensity t. 
Thus the result is in units of the mean waiting time in 
M/M/1 queue. Here the correlation parameter p is the 


parameter. Fixed 8 = 0.5 
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4b. Smoothed, estimated mean waiting time in the correlated 
MDXMD/1 queue divided by the known mean waiting time in 
the M/M/1 queue as a function of the correlation parameter o. 
The result is in units of the mean waiting time in M/M/1 
queue. Here the traffic intensity t is the parameter. 
Fixed 8 = 0.5. 
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D. ESTIMATED DISTRIBUTION OF WAITING TIME W 
The simulation program was set up to enable us to 
look at statistics on W, and We across the replications 
at various n's. In addition the file system enabled us 
to look at joing properties of the samples at two different 


n's, and in particular at correlations. The reason for 
this is that a sample of size 500 is rather small when one 
wants to look at the steady state distribution of Wie 
namely W. One could increase the number of replications, 
but this 1S constrained by the size of the computer. 
However, Since the simulation has to be carried out well 
beyond the set in of the steady state for safety sake, 

one might as well pool samples of waiting times along 

the sample path if they are sufficiently far apart to be 
uncorrelated, or approximately so. Then instead of looking 
at a sample size of 500, we may pool waiting times at 10 


steps of n, Ny, Nor -- to make a sample size of 5000. 


alae 
Siemeche waiting time distribution for a particular run can 
be shown as in fig. 5a, 5b; these are actually histograms 
Showing many of the estimated sample statistics which 
might be of interest. 

In each step of the simulation runs which were performed 


the correlation between waiting times 1000 arrivals apart 


were obtained and found to be negligible, (see Table 3a, 3b), 


a2 
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Table 3b. Correlation Matrix of Positive Waiting Times 


(No Zeros) 


tomavltreorene lnatces nf 12m the 


Correlated Queue. 





thus justifying the pooled samples. Of course a more effi- 
cient use of the data would have been to calculate a correlo- 
gram of estimated serial correlations along each (stationary) 
sample path and average over the five hundred replications 
to obtain a more precise estimate of each serial correlation 
together with a variance estimate of that estimate. 

Another reason for looking at the samples We. = Vie 3, 
500 at different values of n is to assess the convergence to 
a steady state in terms of the whole distribution of ec 
rather than just its mean value. There are a number of ways 
of doing this display and a particularly simple and graphic 
way 1S to give box-plots for the successive samples on the 
Same scale. Unfortunately, this proved to be difficult with 
presently available software;; thus in Figs. 6a-g we give 
box plots of Ww, '5) for n = 1,000(1,000)70,000. The scales 
in these figures are unfortunately not commensurate, but the 
convergence to steady state can be seen. Note that unlike 
the Wi(i)'s, the variance of the wiii's are not decreasing 
with n; in fact that should converge to var(W). If some of 
the distributional detail in the plots is overwhelming, one 
can look solely at the *'s in‘the boxes, which are the median 
values in the successive samples. Since these samples are, 
as we have seen, approximately uncorrelated, the *'s could 
have been smoothed. 

The plots in Figures 6a-g are tied to zero by the occurrence 
of zero waiting times (on the average (l-t)% of the sample). 


Thus plots of the positive waiting times are given in Figures 7a-g. 
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Figure 5a. Histogram and sample statistics for combined waiting 
data in the simulated, uncorrelated M/M/1 queue, i.e. all 
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Histogram and sample statistics for combined waiting time 


data (with 2,293 zero waiting times removed) in the 
simulated, uncorrelated M/M/1 queue, i.e. all W 
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Figures at top and bottom at left are maximum and 
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Figure 7e. Box plots of the sample waiting times Wy(j), j=1,...,500 
with zero waiting times removed for values of 
n = 47,000 to 50, 000° in steps of 1,000. 
Correlated queue; t = 0.7/5, 8 = Q. * o = 0.7/5. 
Figures at top and bottom at left are maximum and 
minimum sample values. 
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Figure 7f. Box plots of the sample waiting times Wy(j), j=1,...,500 
with zero waiting times removed for values of 
n = 51,000 to 60, ae in steps of 1,000. 
Correlated queue; t = 0.75, 6 = 0. 5. Oa URS: 
Figures at top and a. at left are maximum 
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Another question of interest is how the parameters po and 
t affect the waiting time in the MDxMD/1 queue, and how 
these distributions differ from the exponential distribution 
of the positive values of W, for the uncorrelated M/M/1 
queue. There are several graphical techniques for doing this, 
and some of the results are as follows. 

ee Bos plots are used to compare the distribution of 
waiting times at 8 = 0.5 for several fixed values of traffic 
intensity t as the correlation parameter p is varied. These 
are shown in fig. 8a-d and complement Figure 3b also for 
fixed values of the correlation parameter p, the variation 
of the eo pation of W is shown in Figs. 9a-f. Note that 
in these figures 10 approximately uncorrelated samples in the 
steady state have been combined. 

b) The box miler ane not useful for examining departures 
from exponentiality, and for Beate cose the histograms and 
sample statistics from HIST F were utilized. For the grid 
of values of p and t discussed in ‘V.B) the coefficients of 
variation, skewness and kurtosis were extracted from the 
HISTF outputs (no zeros), and these are compared to the 


exponential values C(k) = 1.0, = 200) > a. 6. 00m tne tabler 4. 


Va 
It is clear that for small po & t the distribution is less 
skewed than the exponential, but as 90 and t increase the 
distribution becomes positively skewed and is certainly non- 


exponential. Thus for p = 0.9 and t = 0.9 we have C(x) = 1.1325, 


ie, = Ze 3289, {oe 8.6449. One explanation for this is that for 
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coefficient of variation - 1.0 


'O 
t 
ll 


P> = skewness - 2.0 


P. = kurtosis - 6.0 


Table 4. Deviations of the Estimated Sample Coefficients 
of Variation, Skewness and Kurtosis for the 
Correlated Queue Waiting Times (Without Zeros) 
from the Known (Exponential) Values for the 


M/M/1 Queue. For various p and t values. 
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Oo small, the service time is decreased by the correlation 
when many people arrive. The distribution is therefore 
tighter than for the M/M/1 case. However as p increases the 
service times become highly correlated through the cross- 
coupling, and this effect, which increases the mean waiting 
time and the skewness of the distribution, is dominant. 

c) Another question of interest is whether for large t 
the heavy traffic approximation holds (Jacobs, 1978) and 
the distribution of W is approximately exponential. This is 
difficult to examine because the time for the queue to reach 
steady state is inordinately high for large t. Thus one case, 
t = 0.995, p = 0.25 was studied and the HISTF outputs are 
Shown in Figures 10a & 10b. The plots and the sample statis- 
tics still show a large amount of underdispersion relative 
to an exponential distribution. Thus convergence to an 
exponential distribution as t > l, 1f it occurs, must be very 
Slow. 

Plots for an extreme case, 9 = 0, 8 = 0, are shown 
for various values of t from 0.25 to 0.99 in Figures lla-e. 

d) The question of judging when the transient in the 
Simulation has died out is very complex and difficult, just like any 
question of detecting a signal in noise. In particular it is 
to be stressed that the more graphic output one has, the better 
off one is. This is another reason for looking at the W, (3) 
samples as well as the Ge samples. 

In Figure l2a We is shown for the MDxMD/1 queue 


with t = 0.99, 9 = 0.25. The judgement was made, and published 
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in Table l, that equilibrium was reached by the time of the 
40,000th customer and this appears reasonable from the 

Bugure, which goes out to n = 280,000. The plot of cm 

is given in Fig. 12b; superposing it on Fig. 12a is a help to 
visualization. It of course has a larger transient than Ww 
Beeause 1t drags in all the previous annual times, so that the 
Brot of Ww gives a better picture of the real steady state. 

Of course, as in any real simulation, E(W) is unknown and 

the asymptote in Fig. 12b is a help in showing the convergence. 


Note too that the normality in Ws decreases with n, but that 


OL es does not. 


'"s 


Note too that since we have determined that Ne 
1,000 arrivals apart are almost uncorrelated, local smoothing 
could be applied to the we plot to obtain a better picture. 
Output from a similar simulation for an M/M/1 queue, 
again with t = 0.99 and 500 replications, is shown in Figs. l2c 
& d. The mean of W is known to be E(W) = 24.75. However the 
plot of Ww is well above this from approximately n = 60,000 
to n = 260,000, and since the correlation between we values 
1s not very extensive, this probably indicates that the 
transient is such that E(w.) Geese sgoulfp MOnoOtonically to 
E(W), as one might expect intuitively or from the simulation 
out to n = 240,000, but actually rises above E(W) and then 
declines. It may even have a damped oscillation. The point 
is that all the data should be examined- formal tests of the 
Similarity of the We samples at say W and W would 


100,000 110,000 


probably have indicated convergence (i.e., similar distributions). 
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Another overlaying the curve of oe on that of v 
is a help usually, as is the horizontal asymptotic in Fig. -1l2d. 
Note that the initial transient in the M/M/1 is longer than 
it is in the MDxMD/1 queue. 

Another question which arises is whether one can 
start the simulation in such a fashion as to reduce the length 
of the transient. An attempt to do this is shown in Figs. 13a-b. 
Here an extreme case, t = 0.995, was taken and the initial 
value We was not taken to be zero, but an exponential random 
variable with mean 29.0. The value 29.0 comes from the pre- 
diction formula IV.1. Roughly speaking there appears to be 
faster convergence. The simulation starting from Ne = Q 
takes a long time to converge and was too extensive to 


be handled on the computer. 


is 
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Figure 8a. Correlated queue. Box plots of a sample of waiting times 


Wns (3). j=l,-.-.-,;500; i=1,...,10; with zero waiting times 


removed, for t = 0.25 and five values of o. For each (t,p) 
pair nj is chosen to be large enough for the queue to be in 
Steady state, and 10 approximately uncorrelated samples at 
Ny 2Nos+--+ MQ are combined. 
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Wn; (i); j=1,...,900; i=1,...,10; with zero waiting times 


removed, for t = 0.50 and five values of 0. For each (t,p) 
pair nj 1s chosen to be large enough for the queue to be in 
Steady state, and 10 approximately uncorrelated samples at 

Ny oMos+++sMrQ are combined. 
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Veiom=leee ss ,cO02 1=1,...,10; with zero waiting times 
removed , for t = 0.75 and five values of 9. For each (t,p) 
pair nj is chosen to be large enough for the queue to be in 
Steady state, and 10 approximately uncorrelated samples at 
Ny 2Nos--+sNyq are combined. 
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Figure 9a. Correlated mae Box plots of a sample of waiting times 
Be O=!, 5900; i=],...,10; with zero waiting times 
aeeren- for a= 0.0 and four values of t. For each (t,o) 
pair n, is chosen to be large enough for the queue to be 
in steady state, and 10 approximately uncorrelated samples 
at Ny) Mos-+-sNyQ are combined. 8 = Q.5. 
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Figure 9b. Correlated queue. Box plots of a sample of waiting times 
Wn; j), jJ=l,...,500; i=1,...,10; with zero waiting times 
removed, for 9 = 0.25 and four values of t. For each (t,09) 
pair nj is chosen to be large enough for the queue to be 
in steady state, and 10 approximately uncorrelated samples 
at Ny Nos--+sNyq are combined. 86 = 0.5. 
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Figure 9d. Correlated queue. Box plots of a sample of waiting times 


We iJ); Jals..-.,500; 1=],...,10; with zero waiting times 
rehoved, for p = 0.75 and four values of t. For each (t.p) 
pair nj is chosen to be large enough for the queue to be 

in steady state, and 10 approximately uncorrelated samples 
at Ny2Nos--+-sMyQ are combined. 8 = Q.5. 
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data in the simulated, correlated MDXMD/1 queue; 
t = 0.995, 8 = 0.5, p = 0.25, 1/E(S) = 4.0, 
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Histogram and sample statistics for combined waiting-time 


j.e. all W(5)'s for j=l,...,500 and n = 91,000(1000)100,000. 
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Figure 11b. Histogram and sample statistics for combined waiting-time 
data (with zero waiting times removed) in the simulated, 


correlated MDxMD/1 queue; t = 0.50, 8 = 0.0, p = 0.0, 
ay eeeeosi-e. dil Wij) s for j = 1,...,500 and 
n = 21,000(1000)30,000. 
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Histogram and sample statistics for combined waiting-time 


data (with zero waiting times removed) in the simulated, 
correlated MDxMD/1 queue; t = 0.75, 8 = 0.0, po = 0.0, 


1/E(S) = 4.0, i.e. all W(j)'s for j = 
n = 21,000(1000)30,000. 
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Histogram and sample statistics for combined waiting-time 
data With zero waiting times removed) in the simulated, 


correlated MDxMD/1 queue; t = 0.90, 8 = 0.0, op = 0.0 
Vi Mee4mOne Ice. all Wn(j)'s for gj = 1,...,500 and 
n = 21,000(1000)30,000. 
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VI, VARIANCE REDUCTION IN THE SIMULATION 


eee ths MULTIPLE CONTROL VARIATES METHOD 

Suppose we want to estimate the mean, Hy of a random 
variable X and suppose V is a zero-mean random variable, 
or can be made so by subtracting from V its known mean. 


Then 
Oy. 1) Y = X- av 


is an unbiased estimator for bee and 


var. 2 ) = = o° ~ a*o* =r COV ey jie, 


2 , 2. 
where O,. is the variance of X, oe Tone nemvarnlance OL V. 


ae , y: ? 
Peocaer for oF to be a minimum given Oo 13, anarcov (|<, vi 


the optimal choice for a is readily found to be 


COV XV | 
2 ‘ , 
O T 
V 


ever. 3) a * 


eme@ethen, substituting back in (V.III), the minimum variance 


is 


2 
2a e : 2 © Cen ie | _ 2 _ 2 
8 = oe a = o, (1 See 
V 


E08 





where Be is the coefficient of linear correlation between 
~ and Ve Thus the variance of Y, aS an eStimator of Woe 
is always less than that of X if |cov(X,V)| > 1. Ina 
Simulation X would be the usual estimator of Wye and Y would 
be some random variable generated in the simulation which is 
correlated with x. 

Generally, COV (X,V) ando* are unknown, and we can uSe in 
its stead an estimate derived from the simulation. This 
idea of reducing the variance of an estimator with a control 


variable is readily extended to a vector V of k control 


V, with covariance matrix Q. Now let 


variates Vie veer Vy 


Y = X-a'V 


where a' = (Oy, ne ay) is a row vector of coefficients, 


and standard results from multivariate analysis show that 


oa = 5c + a'Qa - 2a'R 
y Xx 
Here, i = cov{x,V.]. iMac ope micalechoice £Oor a is 


(assuming, of course that Q is non-singular, for otherwise some 
control variates would be redundant). Then the minimal value 


of oO is given by 
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Ce 2 es 
Ce oe R'Q “R 


Let M be the joint covariance matrix of the control variates, V 


Sige x, 1.e., 


o* R" 
x 

M = e 
R Q 


Then it can be seen that with the optimal 


52 = det M 
Y det Q 





For further details see Beja (1968) and Kleijnen (1974). 


be SELECTING THE CONTROL VARIABLES 

In order to simulate the mean waiting time in the MDxMD/1 
queue, several control variables derived from the known pro- 
perties of the M/M/1 queue and the service and interarrival 
processes for the MDxMD/1l queue were examined and combined 
to be the multiple control variables. Note that in the simula- 
tion the M/M/1 queue and the MDxMD/1 queue are run with common 
exponential variates, so that the service and interarrival 
times in the M/M/1 queue are the variables from which the 
EARMA type service and interarrival times in the MD x MD/1 
queue are generated. It is well known that in a single-server 


FIFO queue the difference between the cumulative interarrival 


OS 





times and cumulative service times is a good control for oo 
The additional control variables were meant to give any 
greater variance reduction. In the simulation run with 

o = 0.25, t = 0.50 and 8 = 0,50, Subroutine CONVAR was used 
to compute the correlation between waiting time a average 
waiting time ie in the MDxMD/1l queue and several quantities 
with known means from uncorrelated queue, i.e., the M/M/1 
waiting time, the M/M/l average waiting time, the number of 
arrivals since the last regeneration point ce = 0), sum 

of service time, etc. The controls for i and We are quite 
different and are examined Separately. 


1. Controls for the Waiting Time We 


The control variables from the uncorrelated queue 
which have been selected as being the most correlated to 
a are the number of arrivals since the last regeneration 


(V the waiting time (V5), the average of the previous 100 


1)? 
interarrival times (V3), the average of the previous 100 
interarrival times (V,). The correlation matrix for these 
control variables with Wo is shown in table 5. 

Since the mean of these four variables are not zero, 


the known expected value must be subtracted from each and 


the estimator is 


where 
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V = (V1 1V51V3/Vq) 


ax = gtr 
where 
Ro cov 1x, v1] 
and 
X = waiting time Wa (correlated queue). 


The values of E(V) used, are as follow: 


-_— ey 
ae — odes 
Eve = oy 

E(V,) = u. 


Note that even though V the sum of the waiting times 


aut 


in the MDx MD/1 queue is an average of correlated random 


* 
See Appendix I. 


\ 


Oy 





variables, the variables are exponential and their average 
is still We: 

A controlled estimate for Ww, was obtained using Vie 
Vi1VarVy and the results are shown in Table 6 and Figure 14. 
The 500 replications makes it pessihle to estimate, as 
in Table 5, the values of the components in a*. In Figure 
14 the second, fourth, ... box plots are for the sample of 
Boo controlled Wi(3" Ss, jroo. S00. There ascclearly 
less variability than in the box plots (first, third, etc.) 
for the uncontrolled values W,{5), jee, aoe DOO. successive 
pairs of box plots are for different values of n. 

Table 6 gives, for n = 61,000 (1,000) 70,000 statistics 
of the controlled and uncontrolled Wd). | eee er eraeees 00) 
samples. The values S(MEAN) are to be compared. Thus for 
n = 70,000, S(MEAN), the estimated standard deviation of Woe 
is 0.06375, which is to be compared to the wane 0.04632 for 
the controlled sample average. Thus the ratio of the standard 
deviations is 0.04632/0.06375 = 0.726, and the variance 
reduction is (0), 726) = 0.528. This is not as much as might 
have been expected from the multiple controls, but it reflects 
the difficulty of getting further control variables which are 
not themselves highly correlated with the previous control 
variables. 


Zee COntrois for the Sample-path Average Wee 





Different control variables are needed for We and 


ne Since the former includes all waiting times up to n and 
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the latter 1s local. MThus We will only be correlated with 
local controls, while We will be highly correlated with 
controls which contain the whole history of the sample path. 
Thus 5 control variables from the M/M/1l queue were selected 


to control We time, as follows: 


Vy: number of zero waiting times up until n in the 
M/M/1 queue; 


V5: cumulative interarrival time in the M/M/1 queue; 
V3: cumulative service time in the correlated queue; 
Vy? cumulative servie time in the M/M/1 queue; 

V 


Mean waiting time, ie from the M/M/1 queue. 


The estimated matrix from a simulation run with 500 
meebtcations at p = 0.25, 8 = 0.50, t = 0.50 are shown in 
Table 7. 

The expected values of V which are used to center V 


are as follows: 


E(V,) =n(l-t) , where n = no. of arrivals; 
n 

E(V,) = ) 12 (0.562 aaa where E(X,) = mean interarrival 
i=l + time of arrivals in 

the correlated queue; 

n 

E(V,) = peE(SC.) = Nig, where E(SC,) = mean service time 
i=l c of arrivals in the 


MD x MD/1 queue; 
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E(V,) = }) E(S,) =n S where E(S.) = mean service 
i=1 time of arrivals 
in the M/M/1 queue; 
ive) = — h = 3 3 
2 Tee Yee where Ug = Mean service ELMe 


The results of the simulation run with 500 replications 
and p = 90.75, t = 0.75, 8 = 0.50 are shown in Table 8 and 
Fig. 15. From Table 6 we see that, for n = 70,000, S (MEAN) 
for oe with and without control is, respectively 0.001128 
and 0.001709 and we recall that S(MEAN) is the estimate of 
the standard deviation of T. These values are of course much 
smaller than the corresponding values for We Further the 
ratio of the eoemnaced standard deviations is 
0.001128/0.001709 = 0.660035 and the variance reduction is 
W690) - = 0.4356, a rather healthy value. The variance 
reduction shows up dramatically in Figure 15, which is the 
analog of Figure 14. It is also clear that the controlled. 


W_'s are symmetric, possibly normally distributed, and that 


they have reached steady state. 


mie 





Control 
Variable 


OoG79 0.4518 0.1450 


0.1679 1.0 0.6657 Oman © 
0.4518 0.6657 ee) 0.0330 
=OeL9> =O. 2600 . 053937 


0.1450 Ae Sy 1.0 





Control Variables 


V.: Number of arrivals since last regeneration in 
M/M/1 queue; 


>: Waiting time of M/M/1l queue; 


Average of 100 previous interarrival times; 


V4: Average of 100 previous service times in MDxMD/1 
queue. 


Table 5. Correlation Matrix of Waiting Time W, in the 
MDxMD/1 Queue and the Control Variables; 
p= te 250, and @ = 0.50. 
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Vit. VeeveLuUs ION 


A. In simulation the queue, the box plots, of the sample 
of waiting times, W_(j), and average waiting time, Wi(i), 
are helpful in exploring the convergence to steady state. 
Furthermore the box plots of the controlled waiting times, 
compared with the uncontrolled waiting times, are useful 
when examining just how much variance reduction is obtained. 


Note that the variance reduction attained for We is equivalent 


to cutting down the sample size by at least a half. 


B. An approximation for the mean waiting time for MDxMD/1 


queue in the case of 8 = 0.50 and t > 1 is found to be equal 
jmom0.5 + a times the mean waiting time of M/M/1 queue. 


Thus the approximate waiting time is 


E(W) = zh ul (0.5 + eee 


) 

fe one multiple control variate technique is useful in re- 
ducing the variance of the waiting time and mean waiting time 
when simulating the correlation queue. This is done by simu- 
lating it in parallel with the uncorrelated queue and then 
using the data from the uncorrelated queue to control the 
correlated queue. We found that the effectiveness of the 
variance reduction scheme 1s better for Woe the sample path 


~ 


mean waiting time, than for the waiting time Wie 


Wk Ys 





APPENDIX I 


Derivation of Mean Time to the Last Regeneration Point 
(Zero Waiting Time) in an M/M/1 Queue 


Let X. = the number of customers served in the 
jth busy period of an M/M/1 queue, so the 
x, are 1.i.d. with some distributions 
P. = P(X. = cj) Bo le ag se E(%,) =U, 


Var (X,) = a, € Sia, 


5 
Let N = ) X, be the number of customers in the first 
1=1 
r busy periods 


and let N, be the number of busy periods of length j in the 


first r busy periods, so that 


Now select a customer at random from the first N, let E. 
be the event that the customers are in a busy period of length 
Jj, and let S be the number of customers preceding the selected 


customer in his busy period. Then 


il Pee) Sy NN. /N 
) ( 5 3 af 
ak ; 
By P(S = KES) =f Ome Oy layers ye j=] 

ye 4 

Hence E(S|E.) = = ) OR Ee ee 5 ahve 
J J k=0 é 
1g =i 
E(S) = ) a 5 N./N 
j=1 ; 


sb ase 
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iene, 2 = 
roo J J 


and, 


Lime Ny eS Wiehe wrOodbLliLtty 1, 


roo 


then 


1 ’ a Pp: 
Lig 12 (@) ey ee — a 
roo j=1 u 
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Since the mean number of customers served in a busy period, 


Z 
nn Ls a and varlance a iS ipemses Cox Ssmith, 1961, 
at ro) 
p. 158), then 
iim 9S) 2 eee (t+ t“)/(1 - t)? _ 1) 
ae 2 ‘l-t 1/(lL - t) 
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PROGRAM LISTINGS 
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